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1. Introduction 

Modern density-functional (DF) theory of non-uniform systems [1] continues to yield 
useful insight into the fundamental nature of the liquid-solid (freezing) transition, 
as well as practical approximations for basic thermodynamic properties of the solid 
phase. Treating the solid as a highly non-uniform system, and taking as input known 
structural and thermodynamic properties of the corresponding uniform system, the 
DF approach provides an approximation for the free energy of a solid of prescribed 
symmetry, as a functional of the one-particle density. From the solid free energy follow 
freezing parameters, such as the densities of coexisting liquid and solid phases and the 
Lindemann ratio, as well as various solid-state properties, including the equation of 
state and relative stabilities of competing structures. 

The Ramakrishnan-Yussouff (RY) [2] perturbative DF theory of freezing, which 
is based on a truncated functional Taylor-series expansion of the non-uniform system 
grand potential about that of the uniform liquid, initiated a number of subsequent 
reformulations. Among these are several non-perturbative theories, which posit a 
thermodynamic or structural "mapping" of the non-uniform system onto the uniform 
liquid. Although differing to varying extents in both philosophy and manner of 
approximation, all formulations of the theory require the same essential input, namely 
the direct correlation function (DCF) of the uniform liquid. For classical systems, the 
latter is trivially related to the static structure factor, which is directly obtainable 
from integral equation theories, simulations, or scattering experiments [3]. Numerous 
applications of DF theory to a wide range of bulk and interfacial systems have clearly 
demonstrated its utility as a practical tool for studying classical non-uniform system 
phenomena [4]. 

In recent years, extensions of the DF approach to quantum systems at zero 
temperature have also been developed and applied to freezing of various ground-state 
systems [5], including 4 He [6, 7], electrons [8, 9, 10] (Wigner crystallization of the Fermi 
one-component plasma), and Bose hard spheres [11]. For finite temperatures, a path- 
integral-based extension of the (classical) RY theory [12, 13] has been proposed and 
applied to freezing of helium. However, a general theory capable of describing freezing 
both in the ground state and at finite temperatures is still not at hand. 

Ground-state quantum DF theory is formally similar to the classical theory, but 
a major difference is that it requires as input the quantum analog of the DCF of the 
uniform system. A significant complication for the quantum theory is the lack of any 
simple relationship between the quantum DCF [defined below by (8) and (9)] and the 
static structure factor. Since the latter is generally much easier to compute, in practice 
it is usually necessary to resort to approximate relationships. Unfortunately, this has 
tended to blur the distinction between the issues of the accuracy of the theory and 
the accuracy of the input to the theory. Recently, powerful quantum Monte Carlo 
(MC) methods have been employed to directly compute quantum DCFs for certain 
systems [14, 15], thereby raising the hope that more definitive tests of DF theory will 
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soon be possible. Nevertheless, it remains worthwhile to investigate the utility of 
quantum DF theory using as input the currently more accessible, if somewhat less 
accurate, liquid-state structural data. 

In previous work [11], we extended one version of DF theory, based on the 
modified weighted-density approximation (MWDA) [16], from classical to ground- 
state quantum systems and demonstrated its application to freezing of the Bose hard- 
sphere (HS) liquid at zero temperature. For the liquid-state input data, we used an 
iterative procedure based on (i) the paired phonon analysis (PPA) [17] to solve the 
Euler-Lagrange equation for the optimum pair pseudopotential, (ii) the hypernetted- 
chain (HNC) integral equation [3] to compute the corresponding ground-state energy 
and static structure factor, and (Hi) the Feynman approximation to connect the static 
structure factor to the quantum DCF. In an attempt to compensate for inaccuracies 
in the input data, we simply scaled the DCF by a constant factor at all wave vectors. 
Despite a favourable comparison of predicted solid ground-state energies and freezing 
parameters with available simulation data, the sensitivity of the results to the scaling 
ansatz remained unclear. 

The main purpose of this paper is to make use of significantly more accurate 
data for the energy and structure of the Bose HS liquid to specifically examine the 
sensitivity of quantum DF theory to the quality of the liquid-state input data. In the 
process, we also examine the qualitative influence of particle statistics by considering 
freezing of both Bose and Fermi HS liquids. Aside from being of fundamental interest 
as idealized model systems in which effects of excluded volume interactions may be 
studied in isolation, HS systems also serve as useful reference systems for perturbation- 
theory descriptions of real quantum systems, such as helium and dense neutron matter. 
In Sec. 2 we begin by reviewing the principles of DF theory and the MWDA in the 
context of ground-state quantum systems. Section 3 concerns the structure of the 
liquid and describes (i) calculation of the static structure factor for Bose hard spheres 
via the PPA and an enhanced HNC integral equation that includes four- and five-body 
elementary diagrams [18], (ii) approximate connections between the static structure 
factor and the quantum DCF, and (Hi) adjustment of the Bose liquid energy for 
Fermi statistics using the cluster expansion method of Wu and Feenberg [19]. In 
contrast to our previous study, no scaling or any other modification of input data 
is now employed. In Sec. 4 we outline our application of DF theory to freezing of 
quantum liquids and present results for hard spheres. Considerable improvement 
in the ground-state solid energies results from use of the more accurate liquid-state 
data, and distinct freezing transitions are obtained for both Bose and Fermi systems. 
Quantum statistics are found to significantly influence the freezing parameters, the 
energetically less stable Fermi liquid crystallizing at lower density, and with weaker 
atomic localization, than the Bose liquid. Quantitative discrepancies between theory 
and simulation are discussed in the context of the sensitivity of the theory to the 
accuracy of the input data. In Sec. 5 we summarize and conclude by suggesting 
the need for further improvement on the side of liquid-state theory. Finally, in 
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the Appendix we examine the exact short-wavelength asymptotic behaviour of the 
quantum DCF and possibly important implications for the accuracy of the Feynman 
approximation. 



2. Density-Functional Theory 

At zero temperature the central quantity of interest in DF theory is the total ground- 
state energy E[p], a unique functional of the (spatially- varying) density p(r) that is 
minimized by the equilibrium density [20], i.e., 

Sp(r) 

In general, E[p] may be separated, according to 

E[p] = E id [p)+E c [p) + E ext [p], (2) 

into an ideal-gas energy E id [p] (energy of the non-uniform system in the absence of 
interactions), a correlation energy E c [p] due to internal interactions and exchange, and 
an external energy 

E ext [p\ = J drp(r)(f) ext (r) (3) 

due to interaction with an external potential (f> e xt(f)- One advantage of this separation 
is that for bosons E id [p] is given by the exact expression [11] 

E M = ^jdr\V^pJr)\\ (4) 

where m is the particle mass. Another advantage is that E c [p], although not known 
exactly for non-uniform systems, is amenable to approximation by any of the standard 
DF approximations applied to the classical excess free energy. Here, as in previous 
work [11], we approximate E c [p] by an extension to ground-state quantum systems of 
the modified weighted-density approximation [16]. 

The MWDA is based on the general assumption that the non-uniform system may 
be mapped onto a suitably chosen uniform system [21]. In the context of ground- 
state quantum systems, this implies that the average correlation energy per particle 
of the non-uniform system is equated to that of a uniform liquid e, but evaluated at 
an effective density p, i.e., 

£ ww M/fe(p), (5) 

where N is the number of particles. For classical solids, as well as quantum solids 
obeying Bose statistics, the natural choice of uniform system for the mapping is the 
corresponding liquid. In the case of Fermi systems, however, the effect of particle 
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statistics is known to be quite different in the liquid and solid phases. In particular, 
the exchange contribution to the energy is considerably larger in the uniform liquid 
than it is in the solid, in which exchanges amongst site-localized atoms are relatively 
rare. Comparisons of MC results for solid 3 He and 4 He [22,23], for example, show 
a negligible effect of statistics on ground-state energies. Furthermore, the energies 
for solid 3 He as calculated by Ceperley et al [23], using a MC method to sample 
the square of a properly antisymmetrized wave function, are in good agreement with 
the VMC results of Hansen and Levesque [22], who used an unsymmetrized Jastrow 
trial wavefunction (see Sec. 3). These considerations motivate our assumption that 
within the MWDA the Fermi solid, like its Bose counterpart, is best mapped onto the 
corresponding Bose liquid. This is equivalent to ignoring the symmetry with respect 
to particle exchange of the solid many-body wavefunction, on the grounds that atoms 
in the solid are sufficiently localized about identifiable lattice sites as to be essentially 
distinguishable. The physical interpretation of the MWDA for quantum solids is thus 
the following: The correlation part of the ground-state energy of both the Bose and 
Fermi solids is represented by that of the corresponding Bose liquid of appropriate 
effective density. As this effective density is significantly lower in practice than the 
average solid density (see Sec. 4), the solid correlation energy is in fact modelled by 
that of a relatively low-density liquid, thus properly reflecting the fact that interatomic 
correlations are reduced upon confinement of the atoms to equilibrium lattice sites. 
Clearly in the case of more weakly inhomogeneous quantum systems, such as thin 
films or liquids near interfaces, neglect of exchange effects would be inadequate. An 
interesting open issue for future consideration is then how quantum DF theory might 
be generalized to describe Fermi systems of arbitrary inhomogeneity. 

The effective or weighted density p is now assumed to depend on a weighted average 
over the volume of the system of the physical density, and is defined by 

p=j^ J drp(r) J dr'p(r')w(r - r'; p), (6) 

with p appearing also implicitly as the argument of the weight function w. 
Normalization of w, according to 



/ 



dr'w(r-r';p) = l, (7) 

ensures that the approximation is exact in the limit of a uniform liquid \p(r) — > p]. 
Unique specification of w follows from further requiring that Ef WDA \p\ generate the 
exact "particle-hole interaction" [18] in the uniform limit, i.e., 

/ X2 rpMWDAr ] . 

l im ( ; Nr - [Pl ) =v(\r-r'\;p), (8) 
p(rH P V 8p(r)6p(r') J 11 bl h 1 ' 

where v(\r — r'\;p) is to be interpreted as the quantum analog of the classical Ornstein- 
Zernike two-particle DCF, henceforth referred to as the quantum DCF. Formally, 
equation (8) ensures that a functional Taylor-series expansion of E^ WDA [p] about 
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the density of a uniform reference liquid is exact to second order, and also includes 
approximate terms to all higher orders [16,21]. 

A relationship that proves useful in approximating the quantum DCF (see Sec. 3) 
is obtained by taking two functional derivatives of (2) and using (1), (3), and (4), 
yielding 



where 



-i n /, \ ocp ext (r) 

x (l---l^) = ^FT ( 0) 

is the functional inverse of the static density-density linear response function and 

Y - 1 (\r-r'\-p) = 52Etd[p] (11) 

Xo 11 UP) 8p(r)Sp(r'Y 1 J 

is its non-interacting (free-particle) limit. The Fourier transform of the latter can be 
expressed [using (4)] in the simple form 

In Fourier space, where computations are most conveniently carried out, the weight 
function is given by the analyic relation 

w(k; p) = ^-i-y v(k; p) - 5 kfi pe"(p) , (13) 

which follows directly from (5)-(8). Note that normalization of the weight function 
implies, via (13), the quantum compressibility sum rule 

v (k = 0; p) = 2e'(p) + pe"(p) = mc 2 /p, (14) 

c being the speed of longitudinal sound. This sets a consistency constraint on the 
liquid-state energy and structure. 

In summary, equations (5), (6), (7), and (13) constitute the MWDA for the 
correlation energy of a quantum solid at zero temperature. Practical implementation of 
the theory requires knowledge of the liquid-state correlation energy e and the quantum 
DCF v, the subject of the next section. 



3. Liquid-State Energy and Structure 

For the ground-state Bose HS system, the liquid-state energy and static structure 
factor have been obtained by a method that combines the paired phonon analysis 
(PPA) [17] with integral equation theory. The extraction of the key quantum DCF 
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from the static structure factor is described in detail below. First, however, we briefly 
summarize the PPA. 

The basis of the PPA, as well as of variational Monte Carlo (VMC) methods, is a 
trial ground-state wavefunction of the general Feenberg form 

Mr 1 ,...,r N )=exp[-^2u 2 (r i ,r j )- w 3 (n, Tj, r k ) ], (15) 

i<j i<j<k 

from which the optimum correlation factors u n are determined by minimization of the 
total ground-state energy, according to 



5Ur 



{^b\H\^ b ) 
(iPb\iPb) 



= 0, (16) 



where H is the Hamiltonian. In order to be able to directly compare with results of 
VMC simulations, we mainly concern ourselves here with a trial wavefunction of the 
Jastrow form, which includes only the pair correlation factor, or pseudopotential, 
u 2 (r) = u{r). We have also considered, however, the effect of including triplet 
correlations and briefly comment on this in Sec. 4. The pseudopotential is constrained 
by the general condition that it(r)— >0 as r — > oo, implying vanishing correlations at 
infinite separation, and the condition (specific to hard spheres of diameter a) that 
u(r) — > oo as r — ► (j, ensuring a minimum pair separation of a. In applications to 
solids, VMC methods assume (15) multiplied by a product of one-body (typically 
Gaussian) wavefunctions, each centred on a lattice site. 

By exploiting the formal correspondence between the quantum probability density 
\i(jb\ 2 and the classical Boltzmann factor [24], in which u(r) may be interpreted as a 
classical pair potential, the method uses integral equation theory for classical liquids 
to determine the radial distribution function g(r, [u]) for a given u(r). In previous 
work [11], we used the approximate HNC integral equation, which entirely neglects 
the bridge function. Here, we use a more accurate enhancement of the HNC equation 
(denoted by "HNC+") that approximates the bridge function by four- and five-body 
elementary diagrams. From g(r, [it]), the corresponding correlation energy is computed 
via [17, 24] 

E c {p) = \p J drg(r,[u])[^-V 2 u(r) + ^r)l (17) 

where <f>(r) is the pair potential (for hard spheres in our case). The random phase 
approximation is then used iteratively to adjust the pseudopotential until the energy 
is minimized. Figure 1 illustrates the resulting static structure factor S(k), which is 
related to g(r) by the Fourier transform relation 



S(k) = 1 + p J dr[g{r) — 1] exp(zfc ■ r). 



(18) 



To date there is no known exact relation connecting S(k) to the quantum DCF 
v(k). However, as shown in Sec. 2, v(k) is trivially related to the static response 
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function x(/c). Further, the fact that S(k) and x(/c) may be expressed as frequency 
moments 

P oo 

m p (k)= duj(huj) p S(k,uj), (19) 
Jo 

of the dynamic structure factor S(k,ui) provides a basis for approximations [25]. The 
relevant lowest-order moments (or sum rules) are given explicitly by 

m_i = -^x(k), (20) 
2p 



mo = S(k), (21) 



and 



h 2 k 2 

^ = ~2m~' (22) 
From the rigorous inequality 

/ du S ^ UJ \ l + ahuj) 2 = m-i + 2am + a 2 m t > 0, (23) 
Jo 

variation of the real parameter a to minimize the left-hand side yields the bound [25] 
/7 , 2mS 2 (k) 

m-i(fe) > (24) 

The treatment of (24) as an equality, rather than merely as an upper bound, is the 
so-called Feynman approximation, 



F , 2mS 2 (k) 

m F ^(k) = -^±1. (25) 



From (9), (12), and (20), the corresponding Feynman approximation for v(k) is 

^)-S(^)- 1 )- <26) 

which is plotted in Figure 2 for Bose hard spheres, using the PPA S(k). 

In the long-wavelength (k — > 0) limit, the Feynman approximation evidently 
correctly obeys (14), since S(k) is known to vanish linearly with k, according to [26] 

fik 

However, the PPA - in common with all liquid theories - does not satisfy the 
compressibility sum rule [equation (14)] relating e and v(k = 0). That is, the speed of 
longitudinal sound derived from e differs from that derived from v(k = 0). At k = the 
discrepancy is unimportant, since the MWDA satisfies the sum rule by construction 
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[see (7) and (13)]. At k ^ any remnant of the discrepancy is important only to the 
extent that v(k) may be affected at relevant wave vectors, e.g., the reciprocal lattice 
vectors of the crystal in the application to freezing (Sec. 4). Recently a simple analytic 
modification of v(k) has been introduced [27] to correct the low- A; behaviour. This 
correction, however, would have negligible effect on v(k) at the wave vectors relevant 
to freezing, since the first reciprocal lattice vector of the crystal generally lies at or 
beyond the first minimum in v(k) (see Figure 2). Thus, although in previous work [11] 
the sum rule violation motivated an empirical scaling of v(k), in the present work we 
make no alteration to the PPA input. 

In the opposite short-wavelength (k — > oo) asymptotic limit, it may be shown 
from (26) that v F (k) — > 0, which is identical to the corresponding limit of the classical 
DCF. In the Appendix we consider a recently-proposed extension to the Feynman 
approximation, involving higher-order moments of S(k, u), and use this to analyse the 
exact asymptotic behaviour of v(k). Interestingly, unlike its classical counterpart, the 
exact quantum DCF tends not to zero, but rather to a density-dependent constant. 
The same property has been discussed recently, in another context, by Likos et al [28], 
and may have important implications for the form of the exact quantum DCF in the 
finite-/c range of relevance for the DF theory of quantum solids. 

For the Fermi HS liquid, we have approximated the ground-state energy by 
adjusting the Bose liquid energy for Fermi statistics using the Wu-Feenberg (or 
Pandharipande-Bethe) cluster expansion method [19]. The same approach has been 
adopted previously in VMC studies of liquid 3 He [22] and of the Fermi HS liquid [29] . 
The Wu-Feenberg (WF) method assumes that the Pauli exclusion principle acts only 
as a weak perturbative correction to the strongly repulsive pair potential, and it takes 
as a trial wavefunction ip F the Bose wavefunction tjj B antisymmetrized by a Slater 
determinant: 

. . . , rjv) = iPb(vi, . . . , r N ) ■ det[exp(ifc; • rj)}. (28) 

The corresponding ground-state energy per particle e F may be expressed as a cluster 
expansion about the Bose counterpart e B , according to 

e F = e B + e F + ef + e F , (29) 



where 



e F = le Fl (30) 



ef = 24e F [ dxx\l - \x + \x z ) \S{2k F x) - 1] , (3 1 ) 

Jo 2 2 



and 



s 2=-(-^) e F J dxi J dx 2 J dx 3 x 2 12 S(k F x 12 )[S(k F X23) - l][S(k F x 31 ) - 1], (32) 
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with hp = (3-7r 2 p) 1//3 and ep = (h 2 k F /2m) being the Fermi wave-vector and energy, 
respectively, and with the integrals in (32) extending over the unit sphere. Evidently 
ef is the one-particle correction arising directly from the Pauli exclusion principle for 
non-interacting fermions, while ef and ef arise from antisymmetrization of two- and 
three-particle clusters, respectively. We have computed the first-order correction ef by 
standard numerical quadrature and the second-order correction ef by a simple Monte 
Carlo integration algorithm that randomly chooses triplets of points within the nine- 
dimensional integration volume. Figure 3 shows the resulting relative magnitudes of 
successive terms in the cluster expansion. 

At first glance, the apparent rapid convergence of the first few terms in the WF 
expansion would seem to lend some a posteriori justification to the assumption of 
a small perturbation. Serious doubts have been raised, however, concerning the 
general convergence properties of such cluster expansions [23,30,31]. Brandow [30] 
has identified two distinct antisymmetry effects and argued that one of these, the 
so-called "correlation-exclusion" effect, is treated in a non-converging manner, and in 
fact can be properly accounted for only by including diagrams involving exchanges 
amongst any number of particles up to and including all N particles. The MC 
calculations of Ceperley et a/, for a variety of pair-wise-interacting Fermi systems, 
demonstrate that the WF expansion tends to consistently underestimate the energy, 
particularly for soft-core pair potentials and at high densities. Krotscheck [32] has 
shown that conventional cluster expansion methods tend to be unreliable in the case 
of high densities and long-ranged Jastrow functions, and stressed the importance of 
summing to all orders certain properly combined classes of diagrams (Fermi chains). 
The Fermi HNC method, although summing much larger classes of diagrams than the 
WF expansion, also is plagued by uncertain convergency [32]. Despite the somewhat 
doubtful convergence properties of the WF expansion, we nevertheless employ it here, 
pending a more accurate treatment of the Fermi liquid [33], in order to make contact 
with the available VMC simulation data [29] , and to estimate - at least qualitatively 
- the influence of Fermi statistics on the freezing parameters of hard spheres. 

4. Freezing of Quantum Hard Spheres 

The main steps in applying DF theory to bulk solids and freezing transitions are 
parametrization of the solid density, minimization of the total solid energy with 
respect to the parametrized density, and identification of liquid-solid coexistence by 
construction of a Maxwell common tangent. Below we outline the procedure, which 
is described in greater detail elsewhere [11]. 

Assuming a Gaussian density distribution about the lattice sites of a perfect fee 
crystal, the solid density takes the form 

Ps{r) = (-)§5>xp(-«|r - R\ 2 ), (33) 

71 R 
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where the localization parameter a determines the width of the Gaussians centred on 
the lattice sites at positions R. This choice corresponds to including a Gaussian one- 
body factor in the trial wavefunction of (15). We focus here on the fee crystal structure, 
since solid helium is known from experiment to assume a close-packed structure at zero 
temperature and corresponding pressures [34], and since our previous study indicated 
negligible differences between the relative stabilities of fee, hep, and bee structures. 
The Gaussian ansatz has been widely adopted in applications of DF theory to classical 
systems [1,4]. Simulation studies by Young and Alder [35] and by Ohnesorge et al [36] 
have demonstrated the global form of the classical HS crystal density distribution to 
be close to Gaussian, especially near close packing, although near melting the tails 
of the distribution do exhibit significant anisotropic deviations of roughly 10% from 
the Gaussian form [36]. The only comparable study for quantum systems that we 
are aware of is the Green's Function Monte Carlo (GFMC) simulation of Whitlock et 
al [37] for a Lennard- Jones model of 4 He, which found a spherically symmetric density 
distribution about a given lattice site with only small positive deviations from Gaussian 
behaviour in the tail of the distribution. Since freezing properties are known to be 
dominated by short-range repulsive forces, which are similar in HS and Lennard- Jones 
systems, we can expect the Gaussian ansatz also to be a reasonable approximation 
for quantum hard spheres. Caution is warranted, however, in this application to 
a relatively low-density quantum solid to the extent that the solid energy may be 
sensitive to anisotropy and other deviations of the density from the simple Gaussian 
form. A definitive resolution of the issue could be obtained by the technique of free 
minimization [36] , which avoids entirely the need for density parametrization, although 
we have not attempted such an extensive computation. 
From (13) and (33), (6) now takes the form 



where p s is the average solid density and G the magnitude of the fee reciprocal lattice 
vector G. Given the liquid-state functions e(p) and v(k;p) from Sec. 3, this implicit 
relation for p can be easily solved by numerical iteration at fixed a and p s . From p 



It is important to emphasize that here we map both the Bose and the Fermi HS 
solids onto the corresponding Bose liquid, under the assumption that exchange effects 
in the Fermi solid are negligible. As noted in Sec. 2, this assumption is expected to 
be valid in the solid, where overlap of neighbouring density distributions is negligible. 
In this case, the ideal-gas energy E id takes the simple analytic form [11] 



identical to the form usually assumed in VMC simulations [22] . Assuming no external 
potential, (5) and (35) combine finally to give the approximate total energy per particle 




(34) 



the approximate correlation energy E^ WDA then follows immediately from (5). 



Eid N - t— «, 
4 m 



(35) 



E 



•MWDA 



{ot,p a )/N 



3h 2 



a + e(p(a,p s )). 



(36) 



4 m 
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Note that in the case of helium the mass difference between the 3 He and 4 He isotopes 
gives rise to quite different kinetic contributions to the total energy. For the purely 
kinetic HS systems considered here, however, the energies of the Bose and Fermi solids 
(in units of U 2 /2ma 2 ) are identical by assumption. 

Minimization of (36) with respect to the single variational parameter a (at fixed 
p s ) is illustrated in Figure 4, which shows separately the variation of the ideal-gas and 
correlation energies with a. The linear increase of E^, which results from the kinetic 
energy cost of increasing the curvature of the one-particle wavefunction, strongly 
opposes localization of the atoms about their lattice sites. In contrast, the rapid 
decrease of E c , which arises from a weakening of interactions with reduced overlap 
of neighbouring density distributions, strongly favours localization. At sufficiently 
high density, the competition between E id and E c results in a minimum in the total 
energy at non-zero a, signalling mechanical stabilization of the solid. By varying p s 
and repeating the above minimization procedure, the solid equation of state (E vs 
p s ) is obtained. Figure 5 displays the corresponding weighted density, confirming the 
characteristic trait of the MWDA, referred to in Sec. 2, that p is always lower than the 
average solid density. It should be noted that Figure 5 differs from the corresponding 
figure in [11] (second reference) due to a plotting error in the latter, and that the 
roughness of the curve arises from linear interpolation between grid points of the 
liquid-state v(k). 

Thermodynamic stability of the solid is assessed by comparing the liquid and 
solid energies. Coexistence between liquid and solid, characterized by equality of the 
pressures and of the chemical potentials in the two phases, may be established by 
constructing a common tangent - if one exists - to the energy vs. density curves, the 
coexistence densities occurring at the points of common tangency. The corresponding 
value of a determines the Lindemann ratio L, defined as the ratio of the rms atomic 
displacement to the nearest-neighbour distance in the solid at coexistence. For the fee 
crystal L — (3/aa 2 ) 1 / 2 , where a is the lattice constant. 

Applying the above procedure, we have computed the energy as a function of 
density for the ground-state quantum HS fee crystal and have attempted a freezing 
analysis for the Bose and Fermi systems. In order to assess the sensitivity of the 
theory to the accuracy of the liquid-state input data, we have performed calculations 
for the Bose system using data generated by the PPA with both the simple HNC and 
the enhanced HNC+ approximations restricted to pair correlations (see Sec. 3). The 
PPA static structure factor - and hence the quantum DCF [from (26)] - turns out 
to be practically identical regardless of whether the HNC or HNC+ closure is used. 
The energy, however, differs significantly. Figure 6 shows the Bose HS liquid energy 
density vs number density in the two approximations, together with the corresponding 
fee crystal energy density predicted by our DF theory. Also shown, for comparison, 
are the VMC simulation data of Hansen et al [22]. In both cases, the existence of a 
common tangent to the liquid and solid curves confirms the occurrence of a freezing 
transition. The corresponding freezing parameters are given in Table 1, where it 
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is seen that in comparison with simulation the theory considerably overestimates the 
coexistence densities as well as the density change, and underestimates the Lindemann 
ratio. Evidently, use of the more accurate HNC+ approximation significantly improves 
agreement with simulation for both the liquid and solid energies. In fact, the liquid 
energy agrees almost perfectly up to pa 3 = 0.4, well beyond the freezing density. 
Nevertheless, the coexistence densities are slightly more accurate within the HNC 
approximation, highlighting the delicate dependence of the freezing parameters on 
the detailed form of the equation of state. Despite quantitative discrepancies in the 
predicted freezing parameters, it is important to point out that the prediction of a 
freezing transition itself is a significant qualitative result, considering that the second- 
order RY version of DF theory gives no freezing transition for the same input data [6, 7] . 

To test for sensitivity of the results to the assumption of pair correlations in the 
trial wavefunction of (15), we have repeated the calculations using PPA input data 
that include triplet correlations. Although the liquid S(k), and hence the quantum 
DCF, is essentially the same as with pair correlations, the liquid energies are slightly 
lowered. This results in a consistent lowering also of the solid energies and a slight 
reduction in the densities of coexisting liquid and solid, indicating that comparisons 
with highly accurate simulation data, not restricted to pair correlations, should 
properly incorporate triplet correlations. In comparison with more accurate GFMC 
simulation data [38], however, the predicted solid energies and freezing densities are 
still significantly overestimated. For clarity, it is important to distinguish between 
higher-order correlation factors in the trial wavefunction of the liquid-state theory 
and higher-order DCFs in the approximate energy functional in the DF theory of the 
solid. Although here only the pair DCF is explicitly invoked through (8), in principle 
the theory could be extended to include also higher-order DCFs, defined as higher- 
order functional derivatives of the energy functional, as in fact has been achieved 
already, at the level of the triplet DCF, in the classical theory [39]. 

The same procedure as above has been applied to the Fermi HS system, in this 
case though, using only the more accurate PPA-HNC+ input with pair correlations. 
Figure 7 shows the Fermi liquid and fee crystal energy densities together with the 
Bose liquid energy and the corresponding VMC data of Schiff [29] for comparison. 
We reiterate here that the Fermi liquid energy is approximated by the sum of the 
Bose liquid energy and the WF corrections shown in Figure 3. Of course, because the 
VMC data also are based on the WF expansion for the liquid and on the neglect of 
exchange effects in the solid, the comparison, although consistent, is not a test of the 
validity of these two approximations. The predicted Fermi system freezing parameters 
are listed in Table 1, where the relative change of density upon freezing is seen to be 
about 4 %. Unfortunately, the scarce available simulation data [29] do not permit 
an accurate determination of freezing parameters for comparison. For convenience, 
the predicted energy densities for both Bose and Fermi systems from Figures 6 and 7 
are listed numerically in Tables 2 and 3. Note that the higher energy of the Fermi 
liquid compared with that of the Bose liquid simply reflects the destabilizing effect of 
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Fermi statistics (exchange) on the liquid phase and naturally leads to a lowering of 
the liquid and solid coexistence densities, as is observed experimentally in the case of 
3 He and 4 He [34]. Considering the sensitivity of the freezing parameters to details of 
the equation of state, however, quantitative predictions must await a more accurate 
treatment of the Fermi liquid [33]. 

The consistent over-prediction of solid energies points to a systematic error either 
in the theory or in the input to the theory. To investigate this issue, it is instructive 
to compare the present results with those found previously by a similar analysis [11]. 
In the previous study, empirical scaling of the liquid-state DCF by a factor of roughly 
1.2 considerably improved the freezing parameters. Similar conclusions have been 
reached by other workers [6, 7] who applied the RY version of DF theory to 4 He. On 
the assumption that the MWDA is basically reliable, at least for HS systems, this 
observation would suggest that the true v(k) is considerably more structured than 
that predicted by the current PPA-based liquid-state theory. Indeed, calculations and 
comparison with experiments for liquid 4 He have shown [40,41] that the Feynman 
approximation tends to underestimate the structure of the static response function, 
especially in the crucial region near the roton maximum. This is evidence that the 
incorrect asymptotic behaviour of v F (k) (see Appendix) may be adversely affecting 
the relevant finite-fc region. In quantum freezing, the /c-range of practical interest is 
considerably narrower than in classical systems. This is because the quantum solid 
has a much lower density near the transition than the classical solid, and hence a 
much wider Gaussian density distribution. Included in Figure 2 are the positions 
and relative weights of the first few fee reciprocal lattice vectors, illustrating that the 
region of the first minimum, and to a lesser extent the second minimum, of v(k) is 
crucial in determining the weighted density and thus the energy of the solid. The fact 
that S(k) is relatively insensitive to the type of closure approximation used (HNC or 
HNC+) suggests that it is the Feynman approximation [equation (26)], in converting 
S(k) to v(k), that is most likely responsible for any significant error in v(k). Finally, 
it should be noted that any remaining high-density error in the PPA calculation of the 
input liquid energy will have negligible effect on the solid energy, since the MWDA 
always maps the solid onto an effective liquid of lower density (see Figure 5). Thus, 
the only option for improving the theoretical results appears to lie in the input of a 
more accurate liquid-state DCF, respecting the correct asymptotic limit. 

5. Summary and Conclusions 

By way of recapitulation, we have applied density-functional theory of non-uniform 
systems to quantum hard-sphere solids at zero temperature, and thereby studied the 
liquid-solid transition of Bose and Fermi systems. Mapping both the Bose and Fermi 
solids onto the corresponding Bose liquid through the use of the modified weighted- 
density approximation, ignoring exchange effects in the Fermi solid, and taking liquid- 
state input data from an accurate paired phonon analysis calculation coupled with 
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the Feynman approximation, the solid energy and freezing parameters have been 
computed with no adjustable parameters. The qualitative influence of Fermi statistics 
on the freezing transition also has been considered by using the Wu-Feenberg cluster 
expansion method to approximate the effect of antisymmetry on the Fermi liquid 
energy. Compared with the Bose liquid, the energetically less stable Fermi liquid 
freezes at lower density into a lower-density solid with a higher Lindemann ratio. 

The ground-state energies obtained using liquid-state input data of varying 
degrees of accuracy are generally in satisfactory qualitative agreement with available 
simulation data. Furthermore, distinct liquid-solid transitions are always obtained, 
demonstrating a certain robustness of the theory. Quantitative accuracy depends 
heavily, however, on the accuracy of the input data. The freezing parameters, in 
particular, are especially sensitive to the detailed form of the equation of state. As 
shown previously for Bose hard spheres [11], empirical scaling of the Feynman DCF 
(or static response function) by an appropriate structure-enhancing factor results in 
much closer agreement with simulation than is obtained by using the best available 
unsealed liquid-state data. Experience with 4 He further indicates that the Feynman 
approximation tends to underestimate the structure of the liquid. Finally, analysis 
of the short-wavelength asymptotic behaviour of the exact liquid DCF reveals that 
the Feynman DCF tends to the wrong limit. We conclude therefore that a more 
accurate liquid-state DCF - generated either through alternatives to the Feynman 
approximation or by Quantum Monte Carlo simulation - is still called for to achieve 
a decisive quantitative test of the theory in this demanding application. 
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Appendix. Asymptotic Limit of v(k) 

As discussed in Sec. 3, analysis of the lowest-order frequency moments of S(k,u) 
leads directly to the Feynman approximation for v(k), which is seen [from (26) and 
(27)] to have the asymptotic limit v F (k) — >■ as k — > oo. The moment analysis may 
be straightforwardly extended [41] to include the higher-order moments m2 and 7713. 
The latter may be expressed in the form [25,41] 




(Al) 
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where D(k) is the kinetic structure function, defined by 

D(k) = J dr J dr'cos{k{z-z'))V z V z/ p^{r,r'; R, R')\ {r y }={R , R > } , 



and 



m 3 = 



2m 



+ 



m / 2m z 



drg(r)[l - cos(fc • r)](k ■ V) 2 </>(r), 



(A2) 



(A3) 



with (ex) the mean ground-state kinetic energy per particle and <p(r) the interatomic 
potential. Starting from the rigorous inequality 



I 



* du^^-(l + ahu + bh 2 u 2 f > 0, 
Tiu) 



(A4) 



and minimizing the left-hand side with respect to the real parameters a and 6, yields 
the bound 



m-i(k) > 



m^(k) 



l-A(k)/e(k)' 



where 



and 



A(k) = 



e(k) 



m 2 mi 
mi m 



m 3 /miy _ 2 m 2 
mi Vm / m„J 



/A(k). 



(A5) 



(A6) 



(A7) 



Treating (A5) as an equality, and using (9), (12), and (20), leads to the approximation 



v(k) = 



h 2 k 2 

Amp 



1 



S 2 (k) 



1 



<k)Wkj) 



(A8) 



which in principle represents an improvement over the Feynman approximation. 
Indeed, for 4 He the corresponding density-density static response has been explicitly 
computed by quantum MC and shown to significantly improve the comparison with 
experiment [14]. Since the kinetic structure function is not readily available for hard 
spheres, however, we have not attempted to use (A8) in this paper. Nevertheless, it is 
interesting to consider the asymptotic short-wavelength limit. This is straightforward, 
since it is known [25] that 

2 m 

3h 2 

and since for hard spheres m 3 (/c) may be explicitly expressed as [42] 

/h 4 k 4 \ 2 ffi 2 k 2 \z (h 2 k 2 \i\(P 2, \ ( . \ 

/ IP 1 \ / 1 \ i 

+ (mc 2 ( eK ))(p 2 ( k a)--) , (A10) 
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where P is the pressure and 

_ . . 3 /sinx \ 

Vilx) = cosx , 

x z V x / 



- 3 ^ sinX • (All) 



^ . . 5 rl 2 /smx \ sinxi 

M = Aa + A—- m ')-—]- (A12) 

The various moments appearing in (A6) and (A7) are easily shown to behave 
asymptotically as 

m -> 1, (A13) 

™i -> ^A' ( Al4 ) 
2m 

(h 2 k 2 \* 2h 2 k 2 . , . 

( — ) +77— (^), (A15) 

and 



2m / 3 m 



m 3 



+ (^r) (A16) 



Substituting these limits into (A6) and (A7), we obtain 

m/m - §=$^, (M7) 

which, upon further substitution into (A8), yields 

v(k) - (A18) 
3 p 

Thus, the exact u (fc) approaches a finite negative constant as k — * oo, in sharp contrast 
to the asymptotic vanishing of the Feyman approximation, as well as of the classical 
DCF. It is interesting to note that analogous behaviour has been derived for the 
corresponding function defined for the uniform electron liquid [43]. In fact, as has 
been recently suggested [28], it is quite probably a universal feature of quantum DCFs. 
What this curious asymptotic behaviour may imply for the finite-k behaviour of the 
true v(k) of Bose hard spheres remains to be clarified by a decisive quantum MC 
computation. 
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Tables and table captions 

Table 1. Freezing parameters for the ground-state Bose and Fermi hard-sphere 
systems: Liquid and solid (/cc-crystal) coexistence densities, p l and p s , density change 
A/9, Gaussian width parameter a, and Lindemann ratio L. 

p ; CT 3 p s <r 3 Apcr 3 acr 2 L 

Bose Hard Spheres: 

Theory (PPA-HNC) 0.351 0.365 

Theory (PPA-HNC+) 0.359 0.371 

Simulation 0.23±.02 0.25±.02 

Fermi Hard Spheres: 

Theory (PPA-HNC+) 0.313 0.326 0.013 10.24 0.235 

Table 2. Liquid-state numerical data corresponding to Figures 6 and 7: Predicted 
ground-state energy densities over a range of number densities p for the Bose hard- 
sphere liquid (E B /V) in the HNC and HNC+ approximations and for the Fermi hard- 
sphere liquid (E F /V) in the HNC+ approximation. 



pa 3 


E B /V (HNC) 


E B /V (HNC+) 


E F /V (HNC+) 


0.02 


0.0043 


0.0044 


0.0068 


0.04 


0.0205 


0.0204 


0.0290 


0.06 


0.0539 


0.0524 


0.0702 


0.08 


0.1099 


0.1052 


0.1345 


0.10 


0.1945 


0.1842 


0.2266 


0.12 


0.3141 


0.2954 


0.3521 


0.14 


0.4760 


0.4452 


0.5174 


0.16 


0.6883 


0.6411 


0.7295 


0.18 


0.9601 


0.8909 


0.9964 


0.20 


1.3018 


1.2036 


1.3271 


0.22 


1.7252 


1.5889 


1.7313 


0.24 


2.2434 


2.0574 


2.2198 


0.26 


2.8713 


2.6210 


2.8045 


0.28 


3.6257 


3.2925 


3.4981 


0.30 


4.5252 


4.0857 


4.3147 


0.32 


5.5906 


5.0159 


5.2693 


0.34 


6.8451 


6.0996 


6.3781 


0.36 


8.3141 


7.3546 


7.6586 


0.38 


10.0257 


8.8000 


9.1294 



Table 3. Solid-state numerical data corresponding to Figures 6 and 7: Predicted 
ground-state energy densities over a range of number densities p for the hard-sphere 
solid (E/V) - independent of statistics - in the HNC and HNC+ approximations. 



pa 3 


E/V (HNC) 


E/V (HNC+) 


0.26 


3.3346 


3.0779 


0.28 


4.0648 


3.7279 


0.30 


4.8955 


4.4608 


0.32 


5.8279 


5.2754 


0.34 


6.9521 


6.2449 


0.36 


8.2973 


7.3914 


0.38 


9.8544 


8.7044 


0.40 


11.6209 


10.1785 


0.42 


13.5754 


11.8113 


0.44 


15.7220 


13.5863 



0.014 14.38 0.206 
0.012 13.39 0.214 
0.02 6.25 0.27 
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Figure captions 

Figure 1. Static structure factor S(k) vs. wave vector magnitude k 
for Bose hard-sphere liquids of three different densities, generated by 
the PPA with the HNC+ approximation (see text). 

Figure 2. Quantum direct correlation function v(k), corresponding 
to S(k) in Figure 1, computed from the Feynman approximation 
[equation (26)]. Vertical bars indicate positions of fee reciprocal lattice 
vectors at average solid density p s cr 3 = 0.4, for which the Gaussian 
width parameter aa 2 = 15.8. Bar heights, arbitrarily normalized to 
50, are proportional to the product of degeneracy and the Gaussian 
factor in (34), and thus represent relative contributions to the weighted 
density p. 

Figure 3. First three terms in the Wu-Feenbcrg cluster expansion, 
representing Fermi-statistics corrections to the Bose hard-sphere liquid 
energy per particle (units of K /ma 2 ). 

Figure 4. Hard-sphere solid energy density (units of h 2 /ma 5 ) vs. 
Gaussian width parameter (solid line) for average solid density p s a 3 — 
0.37. Shown separately are ideal-gas energy (short-dashed line) and 
correlation energy (long-dashed line). 

Figure 5. Weighted density vs. average solid density at the energy 
minimum for the hard-sphere fee crystal. 

Figure 6. Bose hard-sphere liquid and solid energy densities vs. 
density, generated for the liquid by the PPA and for the solid 
by density- functional theory (DFT). Dashed lines correspond to 
the simple HNC approximation, solid lines to the enhanced HNC+ 
approximation (see text). Circles and squares denote variational 
Monte Carlo (VMC) data [22] for the liquid and solid, respectively. 

Figure 7. Fermi hard-sphere liquid and solid energy densities vs. 
density, generated for the liquid by the PPA, with the HNC+ 
approximation and Wu-Feenberg corrections (see Figure 3 and text), 
and for the solid by DF theory. Dashed line shows the Bose liquid 
energy density from Figure 6, illustrating the approximate magnitude 
of statistics corrections in the liquid. Circles and squares denote VMC 
data [29] for the liquid and solid, respectively. 
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